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Abstract 


The steady state heat conduction equation representing 
the diurnal variations in the upper atmosphere are 
numerically integrated with extreme ultraviolet radiation 
as the only source of energy input. Two models of vertical 
flow are as summed and horizontal flow is neglected. The 
results yield a maximum at 1800 local time. Other dynamical 
effects are discussed. 
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INTRODUCTION 

The diurnal variation of the upper atmosphere at heights above 200 km 
has been well established by satellite drag observations. During the morning, the 
density increases until a maximum is reached at about 14:00 hours local time. 

Then it decreases rather rapidly in the afternoon and evening, followed by a 
slower rate of decrease during the night till about 04:00 hours local time. This 
is observed for all levels of solar activity. 

Models of this variation have been constructed by various authors (Jacchia 
(1965), Harris and Priester (1962)) relating temperature variations and density 
variations (through the hydrostatic law). In all these models, the temperature 
and density variations are in phase above 200 km. 

The model of Harris and Priester (1962) was a time dependent one which produced 
the observed variations by considerations of the effects heat conduction, solar 
extreme ultraviolet heating and an adjustable artificial heat source. The time 
dependent heat conduction equation in one dimension was integrated assuming 
hydrostatic equilibrium at every instant. Their theory also included a model of 
vertical convective heat transport. Horizontal conduction and convective trans- 
port was assumed to be negligible. It was found necessary to introduce a second 
artificial heat source in order to obtain a maximum at 14:00 hours local time. 

With E. U. V. heating alone the maximum would appear at about 18:00 hours local 
time with an amplitude too large (with the ratio of the diurnal maximum to 



minimum in temperature greater than two, instead of the model produced ones 
of 1.3 or 1.5). Lagos and Mahoney (1967) has obtained similar results with 
E. U. V. heating (though with a smaller variation). They also calculated lati- 
tudinal variations, by allowing for zenith angle variations with latitude, but still 
only integrating the one dimension time varying heat conduction equation. 

In this paper, the three dimensional heat conduction equation is integrated. 
Essentially similiar results are obtained for the latitudinal and longitudinal 
variations as those derived by integrating the time dependent one dimensional 
equation without the extra heat source. It is concluded that dynamical effects 
are sizeable, so that their inclusion will be necessary (and hopefully sufficient) 
to bring about agreement with observation. 


Steady State Heat Conduction Equation 


The time dependent heat conduction equation is 






where is the heat sources and sinks, ttthe conveotive flow velocity, the 
specific heat at constant pressure, K the thermal conductivity and £ the density. 
The diurnal variation appears to be a steady state phenonemon, though its 
amplitude and mean value may vary due to varying levels of solar activity. Thus 
in a frame of reference fixed with respect to the sun the atmosphere appears to 
have a bulge displaced about thirty degrees from the earth-sun line towards the 
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evening side, though the atmosphere may rotate rigidly with the earth. Thus in 
the steady state, we may take the term on the left side as zero, and we are left 
with the three dimensional heat conduction equation. But to an observer on the 
earth, there appears to be dynamical variations. To obtain the appropriate form 
of the heat conduction equation for an observer on the earth, we must first trans- 
form ' equation (1) to a frame rotating with the earth and then take the limit to 
steady state. Let ft be the longitude in the fixed frame and ^'the longitude in 
the rotating frame, and let t' be the time in the rotating frame. Then the 
transformation equations are 

(f > ' ■» - cO t 

t' - t 


where is the angular velocity of the earth. 


The partial derivative with respect to longitude becomes 




4V j 


and the partial derivative with respect to time is 




In the steady state we can let the partial derivative with respect to time go to 
zero and thus obtain the following relationship between the local time variations 
and the longitudinal variations for an observer on the earth 
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(vector quantities may slmiliary be transformed (see Eshenazi (1967)). 

Once one has obtained a solution as a function of t\ one can obtain the 
general solution by replacing the argument t* with Care must be 

exercised when this relationship is used; in particular, one must not try to 
describe transient phenomena using equations based upon this relationship. 
Boundary conditions imposed upon the equations must be consistent with the 
assumption of steady state if such a relationship is used (see Volland (1966). 
The heat conduction equation, upon the assumption of steady state, becomes an 
elliptic type with the periodicity in the longitudinal coordinate (or local time) as 
an essential boundary condition. The solution will be completely determined 
by imposing boundary conditions at the lower and upper boundaries in altitude, 
together with the knowledge of the heat sources and sinks. 


It is to be noted that the solutions of the tidal equations (Chapman and 
Lindzen (1970)) that would correspond to steady state are those for which the 
parameters f and s in tidal theory have the ratio one half. 


For the heating due to the extreme ultra-violet radiation, we take 




( 2 ) 


where & is the optical depth. 



tc V, Ojt) - I \ <r. tf) dK 

As we are interested in along latitudinal variations, the plane earth approxima- 
tion used by Harris and Priester is no longer valid. Care must be exercised 
here in the calculation of the optical depth so that one obtains the correct phase 
of the heating source as a function of latitude and local time. After numerous 
tests of various methods, the following approach was adopted for the calculations 
of the optical depth. Consideration had to he given to the speed of calculation 
and memory size available in the computer. The technique yielded an accuracy 
of about 2% in the calculation of the optical path for critical values of the 
optical path (values between 3 and .002. This gives yielded a maximum numerical 
error of about 7 % In the calculation of the heating. Variations of the density 
along the optical path due to local time variations could be ignored (test showed 
that this caused an error of at most a fraction of a percent. Thus densities 
corresponding to the local time and latitude at which the heating was being calcu- 
lated were used in the calculation of the optical path. The technique is at follow: 

In the numerical integration of the heat conduction equation, altitude steps of 
one km were used. For the calculation of the optical depth, a seven point Gauss - 
Laquerre type integration formula was used, using as a sealing factor the 
seals height of the constituent. Linear interpolation was performed for the 
densities when the aboissas of the integration formula fell between the altitude 
steps. As the scale heights ranged in values from 25 km and up no lost in 
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accuracy was obtained by this procedure. For zenith angles greater than 88 


degrees, the procedure was modified by breaking the integral into two parts. 


The first part covering a range of altitudes from the point at which the absorption 
is being calculated) to a height of one scale height was calculated by a Gauss - 


Chebyshev type integration formula, and then the remainder by 


Gauss -Laquerre type. 


The cooling due to the reradiation on the infra-red from atomic oxygen was 


included. It is sufficient for this type of calculation to take the entire extreme 


ultraviolet flux as composed of two flux lines (see Mahoney (1966), Harris and 


••// 


Priester (1962)) with the following cross sections; 1.5 x 10 cm for nitrogen 


and molecular oxygen, 1.2 x 10 17 cm 1 for atomic oxygen. The total flux used 


in the EUV range was 1.2 ergs /cm* sec and 0.5 erge/cm x sec of ratiation in 


• i 


the Shuns an-Runge region. The initial number densities at 120 km were 4 x 10 / 


cm*, 7.5 x 10 ,C /cm* and 7.6 x i0 /C /cm S for nitrogen, molecular oxygen and 


atomic oxygen respectively* The calculation of the thermal conduction coefficients 


and specific heats was as performed by Harris and Priester (1962). 


For a complete solution of the problem one should also include the Navier- 


Stokes equation together with the continuity equation and the equation of state. 


This will not be attempted here. Instead, we shall use a phenomological approach, 


making reasonable physical assumptions and see what conclusions can follow. 


As the approach to hydrostatic equilibrium is rather rapid compared to the time 
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scale of the phenomenon we shall be interested In (see Thomas (1969)), we shall 
assume hydrostatic equilibrium in the steady state. We are interested In the 
steady state solution so we shall follow the procedure dismissed by replacing 
time derivatives with longitudinal derivatives or vice versa, in a system rotating 
with the earth. 

In this work, we shall not deal with horizontal convective transport of heat. 
For the vertical convective transport of heat, we may consider several modek 
In the construction of these models, one should take care that the interpretation 
in the rotating frame and fixed frame are consistent. The first one we shall 
oonsider is the following. In an earlier paper (Harris and Priester (196$)), it 
was assumed that the diurnal motion of the atmosphere was such that any parcel 
of the atmosphere moved in such a manner that the total pressure above it was 
constant (see Thomas and Chins (1969)). This may be expressed as follows 

t* = O j 

where O /Dt is the total convective derivative, i.e. we have 

> p /it + (3) 

where tt Is the flow velocity. Ignoring horizontal flow and making use of hydro- 
static equilibrium 

*- 
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w® have (where W is the vertical flow velocity and g the acceleration of gravity 
and £ the density. 

Jp 4t * w e. a , 

or 

w * c RT / 3 )^ f - p j 

whsr. R U gas oonatant. Making use of the hydroatatlo law again, we have 


from 




where H is the pressure scale height and die lower boundary and r the 
altitude we are at, we finally have 


- -H Vif( S'**) . 


Thus we have the final form for the heat conduction equation 


Q C P 


( 3 ? jt +W *V*e ) "V- fr %■) 


)ra> St*# 4 1 ^ * c • \ * 


Here £p is the specific heat at constant pressure, Y the vertical distance troa 






time or longitude, the angular velocity of the earth. The calculations of 
the heat sources and sinks have been discussed. The construction of the finite 
difference equation and other numerical procedures are given in the appendix. 

An alternate model of the vertical flow can be derived from the following 
considerations. It has been convenient in studies of fluid flow to assume 


incompressibility and its implication diverenoe-free flow. In the upper 
atmosphere the incompressibility condition cannot be utilized; but it may be, 


because of frictional effects (ion drag, viscosity), that the steady state of the 
fluid can be suoh where the flow may be regarded as divergence-free, i.e. 


V. v 


T7' U * O 
W 
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We shall explore the results of this assumption together with hydrostatic 
equilibrium in the steady state (again ignoring horizontal transport). 


From the continuity equation 


ie /je t v *( ^ ' a 


v*.”* 
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we have 


W 0 - “ ( * 


where We is now the vertical flow velocity. From the perfect gas law we have 


ISlii 


9 




v/. * (-‘S. -*£)/(*■/„- •?«). „„ 

Making use of hydrostatic equilibrium and using the expression forW, we obtain 
(where H is the scale height) 

W 0 - ^ + M /T • (6 b) 

Thus the heat conductivity equation becomes 

Q C V l ^ ( K<3T ^ ^ 3 , (7) 

where Cy is the specific heat at constant volume (as the flow is divergence-free) . 
One may expect different solutions from this form as compared to the earlier form 
as now the local time derivative has a coefficient down by a ratio of the specific 
heats; and thus the solutions should follow the local time variations of the heat 
source more closely than the previous ones. 

The two solutions are compared in Figure 1; where the equatorial exospheric 
temperature is plotted for the same heat input at the equinoxes. It is seen that 
the diverence free solution has a greater amplitude (as expected because of the 
smaller beat capacity), but the maximum is only slightly shifted with respect 
to the solution of the first equation (5). In both cases the maximum appears close 
to sunset not at fourteen hours local time. 

10 



In Figure 2, there is given the solution corresponding to the first case, where 
the convective derivative of the pressure is zero. Curves of constant exospheric 
temperature are plotted for latitude versus local time. It is seen that there is 
deoided latitude variation (as constrasted to the drag results given by Jacchia 
and Slowey (1967)). This is mainly due to the variation of zenith angle as Lagos 
and Mahoney (1967) has found by integrating the time dependent one dimensional 


equation. There appears to be a shift at high latitudes towards noon where the 
heat input is low. 


In Figure 3, the solution is given for the solstice case (in all solutions the 
incident heat flux is the same). Thus these solutions give a decided seasonal 


variation which again is not observed from the satellite drag measurements. 


Discussion of Horizontal Flow : 

The complete solution for the heat transport problem in the upper atmosphere 
would involve the solution of the Navier-Stokes equation for the horizontal 
flow velocities, together with the heat conductivity equation. In addition, any 
coupling with the lower atmosphere has to be considered. Horizontal flow has 
been calculated by several researchers from the diurnal variation as given by 
various atmospheric models. In these calculations gross simplifications of 
the Navier-Stokes equations have been made. The importance some of these 


simplifications can be estimated. 
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The Navier- Stokes equation in the form we shall consider in the rotating 
frame is 


v ~ K & * - - ~ 7 */% * 3 r -~ v ' > ( 8 ) 


where ^ is the density, p the pressure, 3 the acceleration due to gravity, P* the 
ion drag and Rt the viscous drag. 


We shall assume steady state and to make estimates of the order of magnitude 
of various terms, we shall consider the horizontal flow velocity »inng the equator 
during equinox conditions. Furthermore, we shall assume that the vertical flow 
velocity is given by W. 


Letting V be the horizontal flow velocity at the equator the Navier-Stokes 
equation becomes under these assumptions. 


$L 


tL - r X — 
to: it 




dc 


(9) 


At the equator we may ignore Coriolios and the centrifugal term. We shall 
consider what velocities may be obtained from a linearized treatment ignoring 
drag terms, i.e. when we have 


, , ,p, _ _±. s 

t* to.' % 


In order to make an estimate we may replace RT/g by some mean value of the 
scale height R along the equator and integrate to obtain 


4k 
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V-Vo s J/w* H f -w ( } 


where V© , f c are some reference velocity and pressure at some point on the 
equator. As an estimate, we may set In (p/p ) to unity, oil t equal to 4S0 m/sec 
and H to 50 km. Then v - v 0 is about 1000 m/sec, much greater than that given 
by the earth's angular velocity. With such velocities, it is difficult to understand 
how a steady state condition could be obtained. 

If we were to include the non-linear term, i.e. then the solution would 


— , \ 

V . -- r CO (-\ ±M l - V* C l | 


where we would choose the positive square roots in order to obtain the linearized 
result for small pressure gradients. This then implies that for a steady state 


to exist we must have 
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U,r 


Thus to obtain steady state conditions in the upper atmosphere, with the pressure 
gradients that exist, we must include the other terms neglected, i.e. ion drag, 
viscosity etc. 


For the form of the ion drag, we may take 


Fj 
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(see Geissler (1966)), where o: is the ion density which we shall assume to be 
fixed to the field lines, the effective collision frequenoy. Geissler uses a 
value for ^ /# (where N is the number density) of 5.2 x 10 ' cm 4 /sec. 

Stubbe (1968) has calculated an improved value for the case of atomic oxygen 
of about a factor two greater, but we shall use the smaller value in order to 

-5 

underestimate the effect of ion drag. For the ion density, we shall use 2 x io 
cm 5 a rather low value. For an estimate of the velocity, let us choose a value 
given by the earth's angular velocity, for it is only magnitudes comparable to 
this value that can cause appreciable heat transport to effect the diurnal varia- 
tion, i.e. u = . We shall calculate the ratio of the ion drag to the driving force, 

the pressure gradients. The pressure gradient can be approximated by 

/ i ^ 

For we may use a value of 1/4 x l(f , assuming that the logarithm of the 

pressure changes by unity in a half of a day (this is an overestimate). H may be 
taken to be 50 km. Then the ratio of the ion drag to the pressure gradients 


( t ' °)/c '/*r H V p ) 


becomes 4.5/3 > 1. 


It is to be noted that the above ratio is an underestimate. Thus for velocities 
comparable to the earth's rotational velocity, the ion drag force is comparable 
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to the pressure gradients. Thus velocities much smaller than that given by the 
earth's rotational velocity are to be expected. But for the conveotive horizontal 
transport to be effective, it must have a value comparable to the value of the 
looal time derivative in the heat conduction equation. 

Hh ‘ ^ * T 4tr Q4] 


Here is the horizontal velocity, V rt the horizontal gradient of the temperature. 
For the horizontal gradient of the temperature (at the equator) we may take it 
to be (again assuming steady state and for equatorial conditions at the equinox) 
given by the looal time variations and the lefthand side of the above equation 


becomes 


*/«»* 


. * T A 


which as seen implies velocities comparable to the earth's angular velocity in 
order for the above expression to be comparable to the left hand side of (14). 

Another source of energy input into the upper atmosphere are gravity waves. 
It has been suggested that these may propagate upward with sufficient energy that 
it can be comparable to the E. U. V. heating. This implies a dynamical coupling 
between the upper atmosphere and the mesosphere. In the steady state, these 
motions must occur repeatedly. Such an energy input is not inconsistent with the 
assumption of steady state and hydrostatic equilibrium. We can incorporate this 
possibility in the method utilized here. 


Let us again Ignore horizontal flow and assume hydrostatic equilibrium, 
and let U be the vertical upward velocity. Using the following expression 
for NW 



then the total corrective derivative of the pressure may be expressed as 



The heat conductivity equation becomes 



(W-u) +V'k vT- 


where we can regard the terms containg W*U as a correction to equation 
(5). The expression for W identically vanishes at the lower boundary 
(120 Ion) and the value of U at this lower boundary will be determined 
by the propagation of tidal motions from the lower atmosphere to this 
altitude . 

An estimate of this additional term to the right handside of the 
heat conductivity equation can be made as follows. If we assume that 
U does not change appreciably over one scale height, so that we may 
take the difference Vy-U as a constant U* , then we may integrate 


this term over the altitude and obtain an estimate of the peak flux 
due to tidal motions from the lower atmosphere into the upper atmosphere 


Thus wo obtain 


U t , 1 I} d* ^ j 

/ 2 * «* 

where V is the press ire at the lower boundary about equal to 2.6 x 10 ” 2 
dynes /cm 2 . Thus a decimal varying velocity at the lower boundary of 
less than lm/sec would represent an energy input comparable to the E.U.V. 
flux. Harris and Frlester ( 1962 ) have shown that an additional hea£ 
source of this magnitude can be adjusted so as to obtain agreement with 
satellite drag observations. 

Lindzen and Blake (1970) have shown that the semi-decimal component 
of the lower atmospheric tides may propagate to the upper atmosphere 
with a flux in agreement with the above estimate. 

For a fuller description of the steady state of the upper atmosphere 
one needs a model of the tidal motions at the lower boundary together 
with the integration of the Navier -Stokes equation for horizontal flow. 
Though longitudinal flow may be decreased by ion -drag, the numerical 
results presented here (which have large latitudinal temperature gradients) 
indicate that convective transport towards the poles )for which ion drag is 
less important can he significant. But such a flow is strongly coupled to 
the longitudinal flow due to the non-linear and Coriolis terms in the 
Navier -Stokes equation. The calculations are more involved and work toward 
developing a numerical procedure is in progress. In such a procedure, it 
is a prerequisite the successful integration of the three dimensional heat 
conduction equation, e.g. by a method described in this paper. 
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Figure Captions: 


Fig. 1. Equatorial exospheric temperatures for two models of vertical flow. 


Fig. 2. Latitudinal and local time variations of the exospheric temperatures 
for equinox conditions. 


Fig. J. Latitudinal and local time variations of the exorpheric temperatures 
for solstice conditions. 
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APPENDIX 


Derivation of the Difference Equations 


The steady state heat conduction equation in three dimensions in spherical 
coordinates is 

i*- ^’r) + sm 6 7e ( s ( K 

+ <$ = ^ c p( + w 


Here t is the radius vector, © the colatitude, $ the longitude, Q the heat sources 
and sinks as given by equation two, W the vertical flow velocity as given by 
equation (4) or (6b), (> the mass density, K the thermal conductivity. 

The development of the difference equations will start will the procedure 
given in Varga (1962) on the derivation of difference equations for elliptic 
equations in more than one dimension. We divide the three dimensional space 
to increments A\r , , A<f> and denote the various values of Y by the subscript 

j, of 6 by the subscript K, and of ^(or local time) by the subscript or super - 
subsrcipt t. For points not on the boundary let us consider the elementary 
volume which is a six sided paraUapiped centered around the point (Vj , e k , jft 
whose corners pass through the points given by the set of values ( fj + 

A©/2, r-A^/x). The term involving the thermal conductivity may be 
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integrated over the elementary volume and expressed as surface integrals by 
means of Greens theorem. Then the surface integrals are performed by replacing 
the values of the various quantities by their values at the center of the surfaces. 
The term containing is integrated by evaluating its coefficients at the 

center of the elementary volume. All other terms are integrated by assuming 
their values at the center of the elementary volume. Then the first derivitives 
of the temperature with respect to r, e , at are replaced by central 

difference formulae. The values of the physical quantities at }*£ ,**^ are cal- 
culated by taking their mean vlaues at >*. < , Ktrl appropriately. Aside from 
the matrix elements that arise from the condition that the temperature should 

be periodic in longitude the matrix of the coefficients of the difference equation 

♦ 

can be arrange to that one would have tridiagonal matrix whose elements are 
matrices. Also the term should dominate over the term ^(k 

; if the latter term weren't present the heat conduction equation in steady state 

> 

would be of parabolic type. Thus an iterative method in terms of discrete 
approximation to a parabolic differential equation, (Varga (1962)) where wo 
iterate until a periodic solution in longitude is obtained would appropriate. To 
obtain such a method we express the temperature at the point c in terms of 
points at i.e. (writing now the index £ as an superscript) 
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k' 
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and th« longitudinal derivative of the temperature 

i ,T/ *)^ s rT/' 1 ) 


Finally one obtains the following difference equations for points not 


on 


the boundary. 

4^,k f u( k t; + r' T : +/ t *' 
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where 


t t U - -t,v 


and we have taken the thermal conductivity (so as to linearize the equations at 


the point < to be given by 


KCT'/.J = k(tt;, r ) 


Also we have made use of the fact that the major dependence of the thermal 
conductivity on temperature is of the form proportional to the square root 
of the temperature, so we may make use of this to obtain the following approx- 


imation 
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; 

Similary the terms dp are evaluated by using the projected values of 

i 

the temperatures Thus only a slight error is introduced by not iterating 
the ab nations in order to obtain the values of l(» and C f . 

For the boundary condition on t we used that at , the maximum 

value of j , dT /)>v = $ , and at J* t, T(r,e)=T, for all latitudes. The upper 

boundary of the elementary volume for the point is taken to pass through 
Kh and integrating as before one obtains 

l £ r o R^ * (j^* \ ~ Ki# Su\£ H * ^ Air 






A -4 




where by we mean expression (A5)ie evaluated with J & 4M. 

Simillary one obtains from the boundary conditions at the lower boundary 

= 0 R t A * " ( a /, J/., '• 

For the boundary conditions on latitude one may take the temperatures at 
the poles to be the average over boundary of the temperatures at the latitude 
immediately adjacent to the poles. Alternatively one may obtain the temperatures 
at the poles by integrating the one dimenaionrl heat conduction in altitude at the 
poles. Both procedures were used. Both procedures yielded temperatures at 
the poles within ten percent of each other and less than a couple of peroent at 
all other latitudes. The term representing vertical flow was evalutated explic- 
ltily in term of T£ mi and added to Calculation, were performed 
with and without this term and at most fifty degrees difference in the exospheric 
temperatures where obtained. There was no shift in the peak with local time 
of the exospheric temperatures * 

Care had to excerised in choosing the proper increment in longitude and 
latitude, otherwise numerical instability would developed, for the latitude 


nearest the pole. This instability can attributed to the tarn proportional to 
* 

hf in which would change sign unexpectedly, due to the low value of 4 . F 

It was found that this method of deriving difference equations tended to yield f 

L. 

difference equation which were more stable against this type of instability than 

l 

straight forward differencing would yield. 


A-5 


The matrix of the coefficients of equation A1 (including their values at the 
boundary points) can be structured Into block form In which the sublockc are 
labeled by the index J , and the rows and columns within the blocks by the index 
K. Then the matrix considered in block form is tridiagonal and direct inversion 
was performed to obtain a solution (Varga (1962)). 

Calculations were performed at nine and eleven latitude increments with 
an increment in local time of fifteen and thirty minutes. The calculations with 
nine and eleven latitudes agreed well with each other. The resulte for fifteen 
minutes local time increment raise the peak exospheric temperature slightly 
and slightly narrowed the peak as expected. The results presented are for 
fifteen minutes local time and eleven latitudes. The increment in altitude 
chosen was one kilometer. 

The initial temperatures profiles required for the iterations procedures, 

T 4.-< 

and » were set equal to the lower boundary (120 km) value of 335°k 
for all values of j and*. Rapid convergence was obtained, i.e.T k « -T*. 
decreasing by roughly y ki where 96 corresponds to the number of increments 
of longtitude in twentyfour hours, when a fifteen minute integration step was 
used in local time. To speed up the calculation one can used one hour integration 
steps for two days and then reduced it to fifteen minutes. Widely different 
initial profiles were used and they all converged to the same results. 




